********************************************************************************
* Load main dataset
*
* Authors: Fowlie and Reguant
* August 2020
********************************************************************************

use "$buildpath/output/regressions_dataset.dta", clear

********************************************************************************
* Run several specifications
********************************************************************************

* variables to populate
set obs 100000
gen yvar 	= ""
gen eivar 	= ""
gen bacon 	= ""
gen wt 		= ""
gen iv 		= ""
gen model 	= ""
gen control = ""
gen othercontrol= ""
gen deflator= ""
gen naics_fe 	= ""
gen beta 	= .
gen beta_sd = .
gen beta_sd_cl = .
gen beta2 	= .
gen beta2_sd = .
gen beta2_sd_cl = .
gen obs = .

* count	
local i = 1
foreach deflate in "val" "qty" {
	if "`deflate'"=="qty" {
		replace ln_val_prod = ln_qty_prod
		replace ln_val_exp = ln_qty_exp
		replace ln_val_imp = ln_qty_imp
	}
foreach y in "ln_val_prod" "ln_val_imp" "ln_val_exp" {
foreach eivar in "ei" "ei_tot" { 
foreach bacon in "" "bacon" {
foreach wt in "prod_10" "prod_07" "07" "10" "value_07" "value_10" {
foreach iv in "" "iv" {
foreach model in "" "lin_ei" "lin_ei_lag" "lin_ei2" {
foreach control in "" "ln_fep" {
foreach othercontrol in "" "ln_wage" "ln_wage c.te_2007#c.twexb" {
foreach naics_fe in "naics" "naics year" {
	
	di "`i'"

	qui {
	**********************************************
	* Consider comprehensive energy price
	if "`eivar'" == "ei_tot" | "`eivar'" == "ln_ei_tot" {
		replace ln_p_energy = ln_p_energy_ind
		replace ln_fep		= ln_fep_ind
		replace ln_p_energy_ei = ln_p_energy_ind*`eivar'
		replace ln_p_energy_ei2= ln_p_energy_ind*`eivar'*`eivar'
		replace ln_p_energy_te = ln_p_energy_ind*te_2007
	}
	else {
		replace ln_p_energy = ln_p_energy_dir
		replace ln_fep		= ln_fep_dir
		replace ln_p_energy_ei = ln_p_energy_dir*`eivar'
		replace ln_p_energy_ei2= ln_p_energy_dir*`eivar'*`eivar'
		replace ln_p_energy_te = ln_p_energy_dir*te_2007
	}
	replace ln_p_elec_ei = ln_p_elec*`eivar'
	replace ln_p_elec_ei2 = ln_p_elec*`eivar'*`eivar'
	replace ln_p_elec_te = ln_p_elec*te_2007
	replace ln_fep_ei = ln_fep*`eivar'
	
	**********************************************

	**********************************************		
	* Regression	
	
	* command and main independent variables
	if "`iv'"=="iv" {
		local commd 		= "ivreghdfe"
		local price_	 	= "(ln_p_energy = ln_p_elec)"
		local price_lin_ei	= "ln_fep_ei (ln_p_energy_ei = ln_p_elec_ei)"
		local price_lin_ei2 = "ln_fep_ei (ln_p_energy_ei ln_p_energy_ei2 = ln_p_elec_ei ln_p_elec_ei2)"
		local price_lin_ei_lag	= "L1.ln_fep_ei (L1.ln_p_energy_ei = L1.ln_p_elec_ei)"
		local std			= "robust"
		local std_cl		= "cluster(naics)" 		
	}
	else {
		local commd 		= "reghdfe"
		local price_	 	= "ln_p_energy"
		local price_lin_ei	= "ln_fep_ei ln_p_energy_ei"
		local price_lin_ei2 = "ln_fep_ei ln_p_energy_ei ln_p_energy_ei2"
		local price_lin_ei_lag	= "L1.ln_fep_ei L1.ln_p_energy_ei"
		local std			= "vce(robust)"
		local std_cl		= "vce(cluster naics)" 
	}
	local price_model = "`price_`model''"

	* controls
	if (("`control'"=="ln_fep") & "`model'" == "lin_ei_lag") {
		local control = "L1.ln_fep"
	}
	
	* weights
	local weight = "wgt_prod_`wt'"
	if "`y'" == "ln_val_prod"  {
		local weight = "wgt_prod_`wt'"
	}
	else if "`y'" == "ln_val_imp"  {
		local weight = "wgt_imp_`wt'"
	}
	else if "`y'" == "ln_val_exp" { 
		local weight = "wgt_exp_`wt'"
	}
	if "`wt'"=="prod_07" {
		local weight = "wgt_`wt'"
	}
	else if "`wt'"=="prod_10" {
		local weight = "wgt_`wt'"
	}
	else if "`wt'"=="value_07" {
		local weight = "wgt_`wt'"
	}
	else if "`wt'"=="value_10" {
		local weight = "wgt_`wt'"
	}
	
	* outliers
	cap drop outlier
	if "`bacon'"=="bacon" {
		bacon ln_val_prod ln_val_imp ln_val_exp, gen(outlier) replace
	}
	else {
		gen outlier = 0
	}
	
	* running regression -- two sets of standard errors
	di "`commd' `y' `control' `othercontrol' `price_model'"
	quietly `commd' `y' `control' `othercontrol' `price_model' [aw=`weight'] if outlier==0, absorb(`naics_fe') `std_cl'
	
	if ("`model'"=="lin_ei") {
		local sd_cl = _se[ln_p_energy_ei]	
	}
	else if ("`model'"=="lin_ei2") {
		local sd_cl = _se[ln_p_energy_ei]
		local sd2_cl = _se[ln_p_energy_ei2]
	}
	else if ("`model'"=="lin_ei_lag") {
		local sd_cl = _se[L1.ln_p_energy_ei]	
	}
	else {
		local sd_cl = _se[ln_p_energy]
	}
	
	quietly `commd' `y' `control' `othercontrol' `price_model' [aw=`weight'] if outlier==0, absorb(`naics_fe') `std'

	************************************************
	
	**********************************************		
	* Store results and move to next estimation
	replace yvar = "`y'" if _n == `i'
	replace eivar = "`eivar'" if _n == `i'
	replace bacon = "`bacon'" if _n == `i'
	replace wt = "`wt'" if _n == `i'
	replace iv = "`iv'" if _n == `i'
	replace model = "`model'" if _n == `i'
	replace control = "`control'" if _n == `i'
	replace othercontrol = "`othercontrol'" if _n == `i'
	replace deflator = "`deflate'" if _n == `i'
	replace naics_fe = "`naics_fe'" if _n == `i'
	if ("`model'"=="lin_ei") {
		replace beta = _b[ln_p_energy_ei] if _n == `i'
		replace beta_sd = _se[ln_p_energy_ei] if _n == `i'
		replace beta_sd_cl = `sd_cl' if _n == `i'
	} 
	else if ("`model'"=="lin_ei2") {
		replace beta = _b[ln_p_energy_ei] if _n == `i'
		replace beta_sd = _se[ln_p_energy_ei] if _n == `i'
		replace beta_sd_cl = `sd_cl' if _n == `i'
		replace beta2 = _b[ln_p_energy_ei2] if _n == `i'
		replace beta2_sd = _se[ln_p_energy_ei2] if _n == `i'
		replace beta2_sd_cl = `sd2_cl' if _n == `i'		
	} 
	else if ("`model'"=="lin_ei_lag") {
		replace beta = _b[L1.ln_p_energy_ei] if _n == `i'
		replace beta_sd = _se[L1.ln_p_energy_ei] if _n == `i'
		replace beta_sd_cl = `sd_cl' if _n == `i'
	}
	else {
		replace beta = _b[ln_p_energy] if _n == `i'
		replace beta_sd = _se[ln_p_energy] if _n == `i'
		replace beta_sd_cl = `sd_cl' if _n == `i'
	}
	replace obs = e(N)
	local i = `i' + 1
	**********************************************		
	}
	
}
}
}
}
}
}
}
}
}
}

keep yvar-obs
drop if yvar==""
save "$result_regs/regression_results.dta", replace
